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ABSTRACT 

A  new  algorithm  is  presented  using  a  logarithmic  barrier 
function  decomposition  for  the  solution  of  the  large-scale 
multicommodity  network  flow  problem.  Placing  the  com- 
plicating joint  capacity  constraints  of  the  multicommodity 
network  flow  problem  into  a  logarithmic  barrier  term  of  the 
objective  function  creates  a  nonlinear  mathematical 
program  with  linear  network  flow  constraints.  Using  the 
technique  of  restricted  simplicial  decomposition,  we 
generate  a  sequence  of  extreme  points  by  solving 
independent  pure  network  problems  for  each  commodity  in  a 
linear  subproblem  and  optimize  a  nonlinear  master  problem 
over  the  convex  hull  of  a  fixed  number  of  retained  extreme 
points  and  the  previous  master  problem  solution. 
Computational  results  on  a  network  with  3,300  nodes  and 
10,400  arcs  are  reported  for  four,  ten  and  100  commodities. 
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I .  INTRODUCTION 

Multi commodity  network  flow  problems  emerge  when  several 
distinct  commodities  flow  through  a  common  capacitated 
network  and  share  one  or  more  arcs  that  are  subject  to  joint 
capacity  constraints.  The  objective  is  usually  to  find  the 
minimum  cost  flow  given  demands  and  supplies  of  the  com- 
modities. Problems  of  this  type  are  also  referred  to  as 
"multicommodity  capacitated  transshipment  problems"  and 
"multicommodity  flow  problems"  (MCFPs) . 

The  problem  of  optimizing  MCFPs  arises  frequently  in 
logistic  systems.  As  long  as  only  a  single  commodity  is 
involved,  even  large-scale  problems  can  now  be  solved 
routinely  by  specialized  network  codes  that  exploit  the  pure 
network  structure  of  the  problem  (e.g.,  GNET  by  Bradley, 
Brown,  and  Graves  [Ref.  1])  .  Such  solvers  are  not  directly 
applicable  to  the  general  MCFP,  and  general  linear  program 
solvers  are  usually  inapplicable  as  a  result  of  the  large 
constraint  matrix  encountered  with  typical  MCFPs. 

Because  of  the  importance  of  the  MCFP,  much  effort  has 
been  devoted  to  finding  efficient,  specialized  solution 
techniques.  Earlier  surveys  are  given  by  Kennington  [Ref.  2] 
and  Assad  [Ref.  3].  New  and  effective  decomposition  methods 
were   recently   developed   by   Staniec   [  Ref. 4]   and   show 


encouraging  results.  This  paper  extends  that  research  by 
deriving  and  implementing  a  decomposition  for  the  MCFP  based 
on  a  logarithmic  barrier  function. 

A.   STATEMENT  OF  THE  PROBLEM 

In  order  to  formulate  the  MCFP  as  a  mathematical  program 
the  following  notation  is  used  : 

G    =  {I, J}  is  a  network  with  set  of  nodes  I  and  set 
of  arcs  J. 

P    is  the  set  of  commodities   (products)  flowing  on  G. 

i    e  I  is  a  node  of  G. 

j    e  J  is  an  arc  of  G. 

p    €  P  is  a  commodity  flowing  through  G. 

N    is  an  |I|  x  |J|  node-arc  incidence  matrix  for  each 
product  (N-,  =  N2  =  ...  =  Nipi)- 

N    is  an  |I|*|P|  x  |J|*|P|  matrix  with  matrices  N 
along  the  diagonal,  Os  elsewhere. 

A    is  a  |J|  x  |J|'|PI  matrix  (I,I,...,I). 

c    =  (c-,  ,  .  .  . ,  c  j  p  I  )  is  a  vector  of  arc  costs,  length 

Ul  '  IPI  • 
x    =  (x-j,  .  .  .  rXip  I  )  is  a  vector  of  arc  flows,  length 

Ul  '  IPI  • 
b^   is  a  vector  of  joint  capacities  with  length  I J| . 
b1   is  the  vector  (b^,...,b^)  with  length  |JI*|P|- 
b2   is  the  vector  of  supplies  and  demands  for  each 
commodity  with  length  |J|*|P|. 


Then  the  MCFP  may  be  stated  as  follows  : 

(P)      min     ex  (duals)  (1) 

s.t.    Ax  <  b1  (u1)  (2) 

Nx  =  b2        (u2)  (3) 

0  <  x  <  b1  (u3)  (4) 

For  the  constraint  sets  (3)  and  (4)  the  abbreviated  notation 

x  e  F  will  be  used. 

The  stipulation  that  binds  the  flow  of  several  com- 
modities to  joint  capacities  is  given  by  (2)  and  constitutes 
the  set  of  "complicating  constraints".  Without  those  con- 
straints, the  problem  would  reduce  to  independent,  bounded, 
single-commodity  flow  problems.  The  duals  u-,  corresponding 
to  (2)  are  nonpositive.  The  constraints  (4)  are  redundant, 
but  they  ensure  that  x  is  bounded  when  the  constraints  in  (2) 
are  relaxed. 

B.   SOLUTION  METHODOLOGY 

The  two  basic  approaches  to  solving  the  MCFP  (other  than 
a  standard  primal  linear  programming  method)  can  be  charac- 
terized as  either  decomposition  or  partitioning  techniques. 
The  latter  employ  a  special  basis  factorization  within  a 
simplex  algorithm  in  such  a  way  that  portions  of  each 
generated  basis  maintain  characteristics  of  the  pure  network 
flow  problem.  Those  method  are  not  investigated  here.  For 
further  detail  see  Kennington  and  Helgason  [Ref .  5] . 


Decomposition  methods  solve  MCFP  by  using  a  master 
problem  that  coordinates  the  solution  of  single  network  flow 
subproblems .  These  methods  are  attractive,  since  they  may 
require  the  internal  storage  of  only  one  commodity  at  a  time. 
This  approach  can  further  be  divided  into  resource-directive 
and  price-directive  algorithms.  Both  will  be  stated  here  for 
later  reference. 

The  resource-directive  method  uses  a  master  problem  that 
distributes  improving  capacity  allocations  to  the  individual 
commodity  subproblems.  For  this  purpose  MCFP  may  be  written 
as 

(5) 

(6) 
(7) 
(8) 
0  £  y  £  6 1  (9) 

0  <,   x  <  b1.  (10) 

Any  vector  y  satisfying  constraints  (6)  and  (9)  may  be 
interpreted  as  a  "capacity  allocation"  which  apportions  the 
capacity  of  an  arc  across  the  individual  commodities.  The 
inner  minimization  for  a  fixed  y  amounts  to  a  restriction  of 
(P)  .  If  its  solution  is  feasible  in  (P) ,  it  yields  an  upper 
bound  on  the  optimal  solution. 


<P') 

min 

y 

min   ex 

X 

s.t. 

Ay  =  h± 

X 

-  y  <,   0 
Nx  =  b2 

The  subproblems  for  a  particular  allocation  vector  y  are 
of  the  form 

(SP1 (y) )        V(y)  =  min  ex  (11) 

s.t.        Nx  -  b2  (12) 

0  <  x  <  y  (13) 

which   is   a   set   of   single -commodity   minimum   cost   flow 
problems,  solved  independently  for  each  commodity. 
The  master  problem  than  becomes 

(MP1)  min  V(y)  (14) 

s.t.         Ay  <,  b1  (15) 

0  <>   y  <  b1.  (16) 

The   objective   function   in   (14)   is  piecewise   linear  and 

convex.    In  theory,  the  master  problem  can  be  solved  by 

subgradient  optimization  or  by  a  cutting  plane  algorithm. 

Penalty  and  barrier  function  decomposition  are  examples 
of  price-directive  decomposition.  Before  describing  them,  it 
is  worthwhile  looking  at  the  Lagrangian  relaxation  of  (P)  . 
If  the  joint  capacity  constraints  are  placed  into  the 
objective  function  with  multipliers  u,  <,  0,  the  Lagrangian 
dual  of  (P)  is 

(LR)      max    min  L(Ui,x)  =  ex  +  u^ (b^  -  Ax)        (17) 
u^O   x>0 

st.       x  e  F. 


This  problem  corresponds  to  solving   max   LR(u-i)  where 
LR(u1)  is  defined  by 


u1^0 


(LR(u^))  min  ex  -  u-,  (Ax  -  b-,) 

xeF 

=  min  (c  -  u-.  A)  x  +  u-,  b,  .  (18) 

x€F 

Note  that  the  evaluation  of  LR(u-,)  requires  the  solution  of 

|P|  independent  single  commodity  problems.   Furthermore,  for 

any  fixed  u-^  <  0,  LR(u^)  yields  a  lower  bound  on  (P)  and  this 

bound  will  be  tight  if  u-^    is  optimal  in  (P)  : 

LR(u-^  )   =   ex 

LRfu-^)  is  a  piecewise  linear  and  concave  function  that 

can  be  optimized,  in  theory,  by  subgradient  optimization  (for 

example,  see  Fisher  [Ref.  6],  Goffin  [Ref.  7]  or  Sandi  [Ref. 

8]),  or  by  a  cutting  plane  method  (  see  Kelly  [Ref.  9]). 

Optimization  of   LR(u^)   by   a   cutting  plane  algorithm  is 

essentially  equivalent  to  solving  the  MCFP  by  Dantzig-Wolfe 

decomposition  (see  Staniec,  [Ref.  4]). 

However,  even  if  (18)  is  solved  optimally,  it  is  possible 

~        *      ~       * 
that  LR(u-,)  =  ex   for  u-,  ^  u,  .   Furthermore,  we  may  not  be 

able   to   obtain  a   corresponding  primal   solution  that   is 

feasible  to  problem  (P) . 

Solving   LR(u^)   for   any   u-^   generally   allows   some 

constraints  to  be  violated  with  penalty  -u^ (Ax  -  b^).   Other 

penalty  functions  are  possible  which  will,  at  least  in  a 

limiting   sense,   yield   optimal   primal   solutions.     The 


objective  function  in  (P)  can  be  transformed  into  a  nonlinear 
auxiliary  function  Q(h,x)  that  typically  includes  a  polynom- 
ial penalty  term  of  the  form 

(PP(h)  )   min   Q(h,x)  =  ex  +  |  |  |  (Ax  -  bx)  +  |  I !.    (19) 
xgF 

=  ex  +  q(h,x)  (20) 

where  II  '  |  |   indicates  the  t    norm,  and  (Ax  -  b-,)   is  the 

vector  max  (0, Ax  -  b-^),  i.e.  the  vector  of  violations  of  the 

constraint  set  (2) .   The  value  of  h  constitutes  a  positive, 

increasing  sequence  of  penalty  parameters.   Furthermore, 

t  >  1  is  required  for  this  method  to  ensure  convergence. 

The  penalty  function  has  some  attractive  properties,  as 

k  k 

given  in  Luenberger  [  Ref.  10]  :  Let  {  h  }  and  {  x  }  be 
sequences  of  penalty  parameters  and  optimal  solutions  to 
(PP(hk)),  respectively,  with  hk+1  >  hk  and  h°  >  0 .   Then 

Q(hk,xk)  <,   Q(hk+1,xk+1),  (21) 

q(hk,xk)  £  q(hk+1,xk+1) ,  and  (22) 

cxk  <,   cxk+1  .  (23) 

The  algorithm  converges  to  the  optimal  solution.   Thus, 

lim    Q(hk,xk)  =  ex   ,  and  xk  — »  x 
hk->oo 

Since  a  feasible  solution  is  usually  not  obtained  until  final 
convergence  occurs,  the  penalty  approach  is  classified  as  an 
exterior  point  algorithm. 


In  order  to  obtain  intermediate  feasible  solutions  a 
suitable  scaling  or  projection  procedure  may  be  used.  If  we 
model  the  MCFP  to  include  additional  bypass  arcs  which 
satisfy  any  undeliverable  demand  at  a  high  cost  from  a 
supersource,  we  can  assume  that  any  allocation  of  capacity 
satisfying   (6)   and  (9)   is  also  feasible  in  F.    One  pos- 

A. 

sibility  to  obtain  such  capacity  allocations  y  from  a  given 
vector  x  €  F  is  given  by  defining   y  =  CA(x)  as 

r   xpjblj  /  (Ax)j  if  (Ax-b^j  >  0 

*PJ  "1  a  a  (24) 

L*PJ  +  (bx  -  Ax)j  /|P|    if  (Ax-b1)j  <  0. 

A 

Then,  a  solution  of  (SP1  (y)  )  as  defined  in  (24)  yields  a 

A 

valid  upper  bound  at  the  current  iterate  x.   This  allocation 

procedure  was  successfully  applied  by  Staniec  [Ref.  4]. 

The  idea  behind  the  general  barrier  function  approach  is 

just  the  opposite  of  the  exterior  point  penalty  methods. 

Starting  with  a  feasible  solution  which  lies  within  the 

interior  of  the   constraint   region,   a  modified  objective 

function  establishes  a  barrier  against  leaving  the  feasible 

region.   For  the  MCFP,  the  auxiliary  function  with  respect  to 

the  joint  capacity  constraints  then  becomes 

(BP(h))     min   B(h,x)  =  ex  +  b(h,x).  (25) 

x€F 


An  ideal  barrier  term  b(x)  would  take  the  value  zero  for 
all  interior  points  and  infinity  at  the  boundary.  A  sequen- 
tial barrier  function  b(h,x)  omits  this  discontinuity  and  may 

take  the  following  forms  : 

-1  (26) 

b(h,x)  =  h  I .  (b1  -  Ax) . 

or  J 

b(h,x)  =  -  h  Zj  ln(b1  -  Ax) -,  (27) 

where  (b-^  -  Ax)  •  denotes  the  j    component  of  the  vector 
(b-L  -  Ax)  . 

The  first  expression  (26)  is  called  the  "inverse  barrier 
function"  and  the  second  term  (27)  the  "logarithmic  barrier 
function"  (see,  for  example,  Fiacco  and  McCormick  [Ref.  11] 
or  Bazaraa  and  Shetty  [Ref.  12])  .  Both  functions  approach 
the  ideal  barrier  function  b(x)  as  h  — >  0 .  Any  barrier 
algorithm  requires  the  existence  of  an  interior  region,  i.e. 
it  does  not  work  for  equality  constraints. 

The  logarithmic  barrier  function  was  first  proposed  by 
Frisch  1955  [Ref.  13] .  It  was  then  derived  together  with  the 
inverse  barrier  function  from  the  Kuhn-Tucker  conditions  for 
optimality  by  Fiacco  and  McCormick  [Ref.  11].  The  logarith- 
mic barrier  function  has  obtained  recent  attention  in  linear 
programming  due  to  its  fundamental  properties.  Megiddo  [Ref. 
14]  investigates  the  properties  of  a  weighted  logarithmic 
barrier  function  for  general  linear  programs  that  places  the 
nonnegativity  constraints  on  x  into  the  barrier  term  while 
requiring  strictly  positive  values  for  x.    This  approach 


leads  to  a  smooth  path  through  the  interior  of  the  constraint 
region  towards  the  optimal  solution  and  theoretically 
validates  its  use  in  linear  programming.  Gill,  Murray, 
Saunders,  Tomlin  and  Wright  [Ref.  15]  established  a  close 
relationship  to  the  projective  linear  programming  algorithms 
initiated  by  Karmarkar  [Ref.  16]  . 

The  transformation  of  the  objective  function  into  a 
penalty  or  barrier  function  creates  a  nonlinear  programming 
problem  which  requires  an  efficient  solution  method  to  make 
such  transformation  profitable.  Staniec  [Ref. 4]  successful- 
ly applied  the  method  of  restricted  simplicial  decomposition 
(RSD)  in  order  to  solve  the  penalty  function  decomposition. 
The  idea  was  developed  by  Hearn,  Lawphongpanich  and  Ventura 
[Ref.  17].  RSD  solves  any  pseudoconvex  optimization  problem 
with  linear  constraints  by  generating  extreme  points  in 
linear  subproblems  while  a  master  problem  optimizes  the 
original  objective  function  over  a  simplex  derived  from  a 
fixed  number  of  retained  extreme  points  plus  the  last 
iterate,  i.e.  the  last  solution  to  the  master  problem.  The 
implementation  of  this  method  for  the  MCFP  will  be  described 
later  in  more  detail  since  it  proves  useful  for  the  barrier 
decomposition  as  well. 

C.   TEST  PROBLEMS 

A  real-world  large-scale  MCFP  should  be  used  in  order  to 
examine  the  efficiency  of  the  proposed  algorithms.   Staniec 


10 


developed  an  appropriate  problem  that  describes  the 
transshipment  of  conventional  ammunition  from  production  and 
storage  locations  to  overseas  debarkation  points  and  theatre- 
of-war  locations  via  capacitated  road,  rail,  sea  and  air 
transportation  links.  The  product  demands  are  time-phased 
and  the  objective  is  to  minimize  the  weighted  deviation  from 
on-time  deliveries.  The  network  contains  backlogging  arcs 
with  graduated  penalties  and  bypass  arcs  that  satisfy 
undeliverable  demand  at  high  cost.  For  more  details  see 
Staniec  [Ref .  4  and  18]  .  The  same  network  is  used  to  test 
and  compare  the  algorithms  for  four,  ten,  and  100  commodity 
problems.  The  underlying  network  contains  3,300  nodes  and 
10,400  arcs  of  which  about  10%  are  subject  to  non-redundant 
capacity  constraints. 
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II.  THE  CONCEPT  OF  BARRIER  FUNCTION  DECOMPOSITION 

The  original  MCFP  (P)  constitutes  a  linear  program.  It 
is  only  the  size  of  the  problem  that  makes  a  primal  simplex 
algorithm  unattractive  and  computationally  expensive. 
Penalty  and  barrier  functions  were  initially  designed  for 
nonlinear  programming  problems  where  they  proved  useful  in 
converting  a  constrained  nonlinear  optimization  problem  into 
an  unconstrained  nonlinear  problem.  For  the  special  case  of 
the  MCFP,  penalty  and  barrier  functions  provide  good 
decomposition  tools. 

The  barrier  function  decomposition  for  the  MCFP  retains 
the  basic  decomposition  idea  of  the  penalty  function 
approach.  The  overlapping  joint  capacity  constraints  in  (2) 
are  placed  into  the  barrier  term  of  the  objective  function 
in  order  to  enable  the  successive  solutions  of  independent 
single  commodity  problems  in  a  sequence  of  subproblems .  The 
selection  of  the  logarithmic  or  inverse  barrier  function  is 
not  arbitrary  but  has  an  interesting  theoretical  derivation. 

A.   THE  DERIVATION  OF  THE  BARRIER  FUNCTION 

Fiacco  and  McCormick  derive  the  barrier  function  from 
the  Kuhn-Tucker  sufficiency  conditions  for  constrained 
minima  [Ref.  11].  A  good  and  detailed  analysis,  including 
implementations  and  numerical  results,  is  further  given  by 
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Wright  [Ref .  19]  .   The  derivation  shows  that  the  use  of  a 

logarithmic  barrier  function  is  not  arbitrary  since  it  has  a 

very  natural  origin. 

For  the  purpose  of  this  analysis  it  is  convenient  to 

restate  the  problem  (P)  in  the  following  form  : 

(P'  '  )  Min   f  (x)  =  ex  (28) 

xeF 

s.t.         g(x)=b1-Ax>0  (29) 

where  the  objective  function  (28)  has  gradient  vf(x)  =  c  and 
each  constraint  in  (29)  has  gradient  vg  •  (x)  =  -  a^,  the 
negative  of  the  j  row  in  A.  We  associate  the  dual  vari- 
ables H-j.  with  the  constraints  in  (29)  and  note  that  in 
comparison  with  the  duals  u^  of  problem  (P)  ,  \L*  now  takes 
the  opposite  sign  :  u-,  =  -  (1,  . 

Following  the  derivation  of  Fiacco  and  McCormick,  we 
assume  that  there  exist  points  in  the  neighborhood  of  the 
optimal  solution  to  (P'')  such  that  strict  inequality  holds 
for  the  constraints  in  (29)  i.e.,  g(x)  >  0.  Furthermore,  we 
allow  a  perturbation  of  magnitude  h  in  the  Kuhn-Tucker 

sufficiency   conditions   for   optimality.     At   some   point 

*    * 
[x  (h)  ,  Ui  (h)  ]   near   the   optimum   (x  ^u-^  )   the   following 

conditions  have  to  hold  for  h  small  and  for  all  j  €  J: 

{b1    -  Ax)  •  >  0     (primal  feasibility)  (30) 

Mlj  (bi  "  Ax) j  =  h  >  °  (31) 

(perturbed  complementary  slackness) 
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M-ij  ^  0  (32) 

(dual  feasiblity) 
c  -   Zj  |ilj  (-a=>)  =  0.  (33) 

Rewriting  (31)  as  (1-,  •  =  h  /  (b-,  -  Ax)  ■  and  substituting  in 

(33)  yields  : 

c  -  hlj  [(b1    -  Axjj]"1  (-a^)  =  0.  (34) 

Using  the  notation  g •  =  (b^  -  Ax)  •   £  0   for  all  j  e  J  , 

equation  (34)  is  of  the  general  form 

Vg.[x(hH 

vB(h,x)  =  vf(x(h))  -  h  Z  •   =! =   0     (35) 

3        gj[x(h)] 

and  simply  means  that  the  gradient  of  the  objective  function 

B(h,x)  =  f(x)  -hi-  ln(b1  -  Ax)  ■  (36) 

vanishes  at  x(h) .  This  is  the  logarithmic  barrier  function! 
Note  that  no  constraint  qualifications  are  necessary  since 
all  constraints  in  the  MCFP  are  linear. 

Fiacco   and   McCormick   show   that   the   second   order 

sufficiency  conditions  are  also  satisfied  for  B(h,x)   at 

* 
[x  (h)  ,  u-,  (h)  ]  near  the  optimum  x   and  prove  the  existence  of 

x(h)  satisfying  these  conditions.   It  is  also  shown  that  the 

Hessian  matrix  is  positive  definite  for  small  h. 

The  same  authors  obtain  the  inverse  barrier  function 

from  a  modification  in  the  derivation  above  that  enforces 

nonnegativity  in  (32)  by  introducing  a  variable  X    such  that 

x   =  M-i-;-   The  logarithmic  barrier  function  seems  to  be  the 

more  fundamental  approach. 
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B.   PROPERTIES  OF  THE  LOGARITHMIC  BARRIER  FUNCTION 

The  basic  properties  of  the  barrier  function  are  again 
given  by  Fiacco  and  McCormick  [Ref.  11]  as  well  as  by  Wright 
[Ref.  19]  and  can  be  stated  as  follows  for  the  MCFP  : 

k  V 

For  a  decreasing  sequence  {h  }  and  associated  minima  {x  } 

the  following  conditions  hold  : 


B(hk,xk)  <  B(hk  1,xk  1)  (37) 

k 
for  sufficiently  small   h   and  bounded  (b-,  -  Ax)  •,  and 

cxk  <,   cxk-1,  (38) 

Z.  ln(b1  -  Ax).k  <  I-  ln(b1  -  Ax)  .k_1,         (39) 

lim    hk  I-  ln(b,  -  Ax) ik  ->  0,  and  (40) 

hk-»0      J  J 

lim    B(hk,xk)  -+   ex*.  (41) 

hk->0 

The   following   small   example   will    illustrate   these 

properties.   The  problem  consists  of  a  nonlinear  objective 

2 
function  and  a  single  constraint  :  Mm  2x    s.t.  x  >  1  and 

* 
has  the  optimal  solution  x   =  1 .   The  barrier  function 

2 
B(h,x)  =  2  x   -  h  ln(x-l)  has  a  closed  form  solution 

x(h)  =  0.5  +  (0.25  +  0.25  h)  °*5  . 

The  solutions  are  shown  in  Figure  1  for  linearly  decreasing 

values  of  h  from  20  to  0.1.   The  change  in  x  as  a  function 

of  h  and  the  corresponding  logarithmic  term  are  depicted  in 

Figure  2 . 
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Figure  1  :  Trajectory  of  a  Barrier  Function  Optimization 
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Figure  2  :  x(h)  and  b(h,x)  for  h  decreasing 
It  is  interesting  to  observe  that  the  path  of  x(h) 
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towards  its  boundary  x  =  1  is  smooth  and  can  be  well 
approximated  by  a  straight  line  until  B(h,x)  reaches  its 
maximum.  From  that  point  on  the  rate  of  change  in  x(h) 
increases  slightly.  Similar  results  are  obtained  for 
different  objective  functions.  Returning  to  the  multiple 
constraint  barrier  function  for  the  MCFP,  the  maximum  of 
B(h,x)  is  obtained  for  a  value  h'  which  satisfies  the  sta- 
tionarity  condition  Z_iln(b-,  -  Ax(h'))  •=  0.  This  can  only 
occur  if  some  components  of  (b-,  -  Ax(h))  are  less  than  or 
equal  to  one,  i.e.,  some  constraints  are  almost  tight. 
Thus,  we  will  not  be  able  to  observe  the  property  in  (37) 
until  the  final  stage  of  the  algorithm  when  relatively  small 
values  of  h  are  attained. 

The  properties  (38)  and  (41)  are  the  most  important  ones 
for  practical  purposes.  Starting  with  a  feasible  (interior) 
point,  the  algorithm  produces  a  nonincreasing  sequence  of 
objective  values  which  converges  to  the  optimal  value  of  (P) 
and  provides  intermediate  feasible  solutions  along  its  path. 
This  path  of  x  as  a  function  of  h  describes  a  trajectory 
that  has  already  been  studied  by  Fiacco  and  McCormick.  [Ref. 
11]  as  well  as  by  Wright  [Ref.  19]  .  The  existence  of  this 
trajectory  could  be  utilized  in  an  extrapolation  technique 
predicting  the  next  iterate  or  even  x  .  Its  implementation 
for  the  large-scale  MCFP  is  not  investigated  here.   However, 
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we  will  be  able  to  get  a  good  feasible  solution  after  only  a 
few  iterations  without  using  an  extrapolation  technique. 

One  final  property  is  the  robustness  of  the  barrier 
function  method:  Mifflin  [Ref.  20]  showed  that  it  is 
sufficient  to  solve  B(h,x)  at  each  iteration  only 
approximately,  within  a  predetermined  tolerance,  while  still 
achieving  convergence  (at  a  lower  rate) . 

C.   COMPUTATIONAL  CONSIDERATIONS 

Any  barrier  function  technique  requires  an  interior 
starting  point.  It  may  not  be  easy  to  find  a  starting  point 
and  the  performance  of  an  algorithm  is  influenced  by  the 
quality  of  this  initial  solution.  We  utilize  the  capacity 
allocation  mechanism  y  =  CA(x)  for  the  generation  of  an 
initial  starting  point. 

The  transformation  of  (P)  into  a  nonlinear  programming 
problem  requires  an  effective  NLP  solution  methodology.  If 
a  line  search  is  part  of  this  method,  any  discrete-step  line 
search  procedure  along  an  improving  direction  may  cause  the 
evaluation  of  B(h,x)  at  one  or  more  infeasible  points.  This 
requires  extra  precautions  during  the  implementation, 
resulting  in  additional  computation  time. 

Ryan  [Ref.  21]  points  out  that  for  small  values  of  h  the 
auxiliary  function  B(h,x)  becomes  very  "steep  valleyed"  and 
the  gradient  vB(h,x)  can  take  large  values  in  a  small 
neighborhood  of  x(h).    Therefore,  a  termination  criterion 
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based  on  the  magnitude  of  the  gradient  alone  may  be  critical 
as  h  becomes  small;  termination  criteria  based  on  the 
difference  in  successive  objective  values  may  lead  to 
premature  termination.  We  will  consider  both  criteria  and 
take  the  risk  of  not  solving  the  problem  optimally  at  each 
step . 

Another   well-known   difficulty   of   barrier   function 
approaches  is  the  ill-conditioned  Hessian  matrix  of  B(h,x) 
for  small  h  (see, for  example  Bazaraa  and  Shetty  [Ref.  12], 
Wright  [Ref.  19]  or  Ryan  [Ref.  21]).   This  problem  emerges 
only  in  the  final  stage  as  h  ->  0.   For  large-scale  program- 
ming purposes  the  achievement  of  a  solutiom  which  defers 
from  the  optimum  only  by  some  small  value  e  (  e-optimality) 
is  normally  sufficient.    Such  a  solution  may  be  obtained 
before  the  ill-conditioning  becomes  bothersome. 

Another  issue  that  significantly  influences  the  perfor- 
mance of  the  algorithm  is  the  choice  of  the  initial  barrier 
parameter  h  and  its  rate  of  decrease.  Lootsma  [Ref.  22] 
showed  that  the  absolute  difference  |  I  x(h)  -  x  |  |  is  of 
the  order  0(h)  for  the  logarithmic  barrier  function.   Ryan 

[Ref.   21]   uses   this   linear   relationship   to   propose   a 

k      1  — k   0 
generating  relation  of  the  form  h   =  10     h   ,  where 

k  =  1,2,3,...,  and  h   is  positive. 

The  suitable  choice  of  the  initial  value  h   is  a  more 

critical  issue,  since  theoretically  h   can  take  any  value 
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greater  than  zero  and  up  to  infinity.  Intuitively  however, 
h  should  depend  on  the  cost  and  constraint  structure  of  the 
problem.  One  possible  choice  would  be  to  interpret  the 
parameter  h  as  a  scaling  factor  between  the  cost  vector  and 
the  set  of  constraints  placed  into  the  objective  function. 
This  leads  to  the  suggestion  that  we  achieve  initial  balance 
at  the  starting  point  x   such  that  ex   =  h  Z-  In (bi -Ax  )  •  . 

Instead  of  analytically  deriving  h  ,  simply  taking  an 
multiple  a  (a  greater  than  one)  of  the  maximum  cost  value  in 
connection  with  a  constant  rate  of  decrease  has  performed 
well  in  test  problems  : 

h°  =  a  *  cmav  ,  hk+1  =  a  *  hk  for  k  >  1  and  0  <  a  <  1. 

Ill  CI  X. 

In  general,  we  recommend  choosing  h  too  high  rather  than 
too  low  as  the  algorithm  will  more  quickly  adjust  to  high 
values  rather  than  low  values. 
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Ill .   THE  LOGARITHMIC  BARRIER  FUNCTION  DECOMPOSITION 
USING  RESTRICTED  SIMPLICIAL  DECOMPOSITION  (RSD) 

The  basic  idea  of  the  barrier  function  decomposition  for 
the  MCFP  is  to  place  the  coupling  constraints  (2)  into  the 
logarithmic  term  of  the  barrier  function  as  derived  in 
Chapter  2.  The  resulting  formulation  (BP (h) )  constitutes  a 
nonlinear  programming  problem  with  linear  network  flow 
constraints  (3)  and  (4).  Using  the  restricted  simplicial 
decomposition  technique  (RSD) ,  we  will  decompose  the  problem 
into  a  nonlinear  master  problem  and  a  set  of  subproblems, 
which  require  only  the  solution  of  |P|  independent  pure 
network  flow  problems .  The  master  problem  has  a  reduced 
search  space  described  by  a  fixed  number  of  retained  extreme 
points,  which  are  generated  in  the  subproblems. 

Lower  and  upper  bounds  on  the  optimal  solution  to  (P) 
can  be  easily  established.  The  analogy  with  the  penalty 
decomposition  is  interesting  and  worth  pursuing.  The 
solution  method  RSD  used  in  the  penalty  function  decomposi- 
tion also  proves  to  be  effective  for  (BP  (h)  )  and  is  descri- 
bed first. 

A.   RESTRICTED  SIMPLICIAL  DECOMPOSITION  (RSD) 

The   basic   difference   between   a   linear   programming 
problem  and  a  nonlinear  programming  problem  with  linear 
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constraints  is  the  fact  that  the  optimal  solution  will 
normally  not  be  an  extreme  point  of  the  constraint  region. 

The  familiar  Frank-Wolfe  algorithm  [Ref.  23]  takes 
advantage  of  the  specialized  constraints  by  generating  an 
extreme  point  solution  of  the  original  problem  in  a  linear 
subproblem  whose  objective  function  is  the  linearization  of 
the  original  objective  function  at  the  current  iterate 
(given  an  initial  solution) .  A  master  problem  provides  a 
new  iterate  via  a  simple  line  search  between  the  previous 
iterate  and  the  new  extreme  point.  The  main  disadvantage  of 
this  decomposition  algorithm  is  its  susceptibility  to  slow 
convergence,  especially  when  the  line  search  direction 
becomes  orthogonal  to  the  gradient  of  the  objective  function 
(e.g.,  see  Wolfe  [Ref.  24]). 

The  method  of  simplicial  decomposition  is  due  to 
Geoffrion  [Ref.  25],  von  Hohenbalken  [Ref.  26]  and  Holloway 
[Ref.  27].  A  nonlinear  master  problem  replaces  the  line 
search  of  the  Frank-Wolfe  method  by  extending  the  optimiza- 
tion to  the  convex  hull  of  all  extreme  points  generated  in 
the  linear  subproblems .  This  solution  method  relies  on  an 
effective  solution  of  the  nonlinear  master  problem. 
Although  the  master  problems  have  only  simple  convexity 
constraints,  the  increasing  number  of  extreme  points  makes 
the  implementation  of  this  technique  unattractive  for  the 
solution  of  large-scale  programming  problems. 
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The  disadvantage  of  simplicial  decomposition  led  to  the 
idea  of  restricted  simplicial  decomposition  (RSD)  as  develo- 
ped by  Hearn,  Lawphongpanich  and  Ventura  [Ref.  17].  RSD 
limits  the  size  of  the  master  problem  by  fixing  the  maximum 
number  r  of  retained  extreme  points .  The  master  problem 
optimizes  over  the  simplex  of  these  extreme  points  and  the 
iterate  obtained  from  the  previous  master  solution.  Any 
newly  generated  extreme  point  replaces  the  old  extreme  point 
with  minimum  weight  in  the  expression  of  the  current  iterate 
as  a  convex  combination  of  the  retained  extreme  points  and 
the  prior  iterate.  After  solving  the  new  master  problem, 
all  extreme  points  with  zero  weight  can  be  discarded. 

If  r  is  set  to  its  minimum  value  r  =  1,  RSD  specializes 
to  the  algorithm  of  Frank-Wolfe.  For  the  maximum  value  of  r 
(the  finite  number  of  extreme  points)  the  method  represents 
simplicial  decomposition.  The  solution  to  the  master 
problem  becomes  harder  as  r  is  increased,  but  is  rewarded  by 
significant  improvements  in  the  convergence  rate  to  the 
optimal  solution  (see  Hearn,  et  al . ,  [Ref.  28]). 

The  decomposition  of  (P'')  is  formed  in  the  following 

lc  t  h 

manner:  Let  X   denote  a  matrix  in  the  k    iteration  of  the 

master  problem  whose  columns  are  a  set  of  r  extreme  points 

from  F,  let  x     be  the  solution  of  the  previous  master 

problem,  and  let   Xk  =  Xk  U  {xk-1}.   The  master  problem  in 
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terms  of  the  weights  w  at  iteration  k  becomes, 
for  a  fixed  h, 

(MP2k)   min  B(h,Xkw)  =  cXkw  -  h  I-  In (bj  -  AXkw) -(42) 
s.t.     1  w   =1  (43) 

w  >  0,  (44) 

A  k  A  k     A  k  A  k 

which  has  solution  w   in  terms  of  w  and  solution  x   =  X  w 

in  terms  of  x. 

The   subsequent   subproblem  optimizes  the   linear  ap- 

Ak 
proximation  of  B(h,x)  at  x   over  F,  which  is  equivalent  to  : 

(SP2k)       min  VB(h,xk)x  .  (45) 

xeF 

It  is  convenient  to  introduce  the  notation  (b-^  -  Ax)fc  to 

represent  the  vector  (b-,  -  Ax)  whose  components  are  taken  to 

the  t    power.  Then,  (SP2  )  can  be  written  as 

(SP2k)       min  (c  -  h[(b1  -  Axk)_1]T  A) x  .       (46) 
x€F 

Using  the  relationship  in  (31),  we  observe  that  we  obtain  an 

A  A  A].    _1 

estimate  \l^    of  the  optimal  dual  variables  as  \L^=   Mb-^-Ax  ) 
at  each  iteration.    Substituting  into  (46)  yields  a  sub- 
problem  as 

(SP2k)       min  (c  +  £xA)x.  (47) 

xeF 

This  subproblem  resembles  a  standard  Dantzig-Wolfe  decom- 
position subproblem  as  it  would  be  used  in  a  dynamic  column 
generation  approach  to  problem   (P)   with  dual  estimates 

A  A 

u-l  =  -  lli.  This  issue  is  not  discussed  here;  see  Staniec 
[Ref .  4] . 
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Due  to  the  block  structure  of  the  constraint  matrix  N, 
(SP2  )  permits  independent  solutions  for  each  commodity. 
However,  each  solution  to  (SP2  )  is  not  necessarily  a 
feasible  flow  in  (P)  and  of  course,  is  not  an  interior 
point.  As  a  starting  point,  the  barrier  function  approach 
requires  an  initial  interior  point,  which  may  be  obtained  as 
follows.  First  we  solve  a  subproblem  with  h  =  0,  and 
yielding  solution  x  .  Then,  using  the  capacity  allocation 
mechanism  y  =  CA(x  ),   we  further  reduce  the  capacity  by 

A  A 

setting  y'  =  b  *  y    for   0  <  b  <  1.   Finally,  we  solve 

A 

subproblem   (SP1 (y' ) )   to   get   an   initial   interior   point 
solution  to  (P) . 

Starting  with  this  solution  in  the  master  problem  we 
simply  limit  our  search  to  that  part  of  the  convex  hull  of 
extreme  points  plus  the  current  iterate  which  provides  an 
interior  point  solution  as  next  iterate.  This  will  be 
discussed  in  connection  with  the  solution  of  the  nonlinear 
master  problem. 

B.   SOLUTION  OF  THE  MASTER  PROBLEM 

The  reduced  gradient  method  is  a  well-known  approach  to 
solving  a  nonlinear  program  with  linear  constraints  and  is 
used  here  for  solving  the  RSD  master  problem.  The  method 
was  first  proposed  by  Wolfe  [Ref.  28]  and  modified  by 
McCormick  [Ref.  29] .  The  basic  idea  is  the  partition  of  the 
variables   x   into  m  basic  variables   xB  and  n  nonbasic 
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variables  xN  as  done  in  the  simplex  algorithm  of  linear 
programming.  This  induces  a  partition  of  the  constraint 
matrix  A  into  parts  B  and  C,  where  B  is  assumed  nonsingular. 
The  NLP  then  takes  the  form 

min  f(x)  =  f(xB,xN)  (48) 

s.t.   Ax  =  BxB  +  CxN  =  b  (4  9) 

xB,xN  >  0.  (50) 

The  variables  xN  are  regarded  as  independent  variables 
whereas  xB  are  dependent  variables  completely  determined  by 
xN  and  equations  (49) .  Consequently,  the  objective  function 
can  be  considered  to  be  a  function  of  xN  only  and  the 
constraints  reduce  to  the  nonnegativity  constraints  on  the 
independent  variables  and  the  limitation  in  their  change 
that  provides  nonnegative  basic  variables.  This  fact  allows 
the   application   of   a  modified   steepest   descent  method 

accounting  only  for  the  nonnegativity  constraints.    The 

k     .  ... 

reduced  gradient  r   at  iteration  k  is  of  dimension  n-m   and 

computed  as  : 

rk  =  (rBk,rNk)  =  vx  f(xk,xk)  -  vx  f(xk,xk)  B"1C.  (51) 

N  B 

To  find  an  improving  direction  d  such  that  vf(x  )d  <  0,  we 
select  at  each  iteration  k 

-  r  •     if  r  •  £  0 
dN-  =  «/  for  x_i  nonbasic     (52) 

x-r  •      if  r  ■  >  0 
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and 


dBi   =   "  (B  1CdN)i         for  xi  basic.      (53) 


A  new  direction  needs  to  be  computed  as  soon  as  a  nonbasic 
variable  attains  its  zero  level.  If  a  basic  variable 
becomes  zero,  the  partition  must  be  modified.  Also,  this 
method  requires  a  nondegenerate  solution  at  each  iteration. 
A  more  detailed  description  can  be  found  in  Luenberger 
[Ref .  10]  or  Bazaraa  and  Shetty  [Ref .  12]  . 

The  use  of  the  reduced  gradient  method  for  the  solution 

]r 

of  the  master  problem  (MP2  )  results  in  some  nice  simplif- 
ications due  to  the  presence  of  only  a  single  convexity  con- 
straint. There  is  only  one  basic  variable  w-  to  be  selected 
and  the  reduced  gradient  becomes 

r-  =  ex.  -  cxi  +  h[(b1  -  AXw)_1]T  (Ax.  -  Axi)     (54) 
for  the  nonbasic  variables  w  •  .    Computing  the  direction 
component  for  the  basic  variable  w^  reduces  to 

d±  -  -  B"1Cdj  =  -  Ijdj  (55) 

and  for  the  nonbasic  variables  w  • 

f   -  r  .    if  r  •  £  0 


dD  - 


(56) 


^  "w-ir-i   if  r  •  >  0 


The  line  search  in  the  direction  d  is  limited  by  the  non- 
negativity  constraints.    Thus  a  maximum  steplength  a'  is 
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computed  as 

a'  =  min  {-wJc/dIc  I  dk  <  0  ,  for  all  variables  w^}  .  (57) 
Any  move  in  the  direction  d  of  size  a  must  also  satisfy 
b1  -  [AX(w  +  ad)  ]  >  0   for  all  0  <  a  <,   a'  (58) 
or  equivalently 

bj_  -  AXw  >  OAXd      for  all  0  <,   a  <,   a'  .     (59) 
Since  b,  -  AXw  >  0  ,  this  holds  for  all  arcs  1  with 
(AXd)-,  <  0.    If  (AXd)  -i  >  0  for  some  arcs  1,  we  perform  a 
ratio  test  to  select 

a"  <      min  {  (b1  -  AXw) 1    /     (AXd) ±     |  (AXd) ±    >  0)  (60) 
and  set 

amax  =  min  {a'  f°-'  '  }  '  (61) 

The  master  problem  solution  is  summarized  as  follows: 

Step  0  :  Set  w  such  that  lw  =1,  wk  >  0 . 

Select  a  basic  variable:  use  the  largest  w^. 
Step  1  :  Compute  the  direction  d  as  determined  in  equations 

(54)  through  (56) .      If  d  =  0,  stop. 
Step  2  :  Solve   min  B[h,X(w+ad)]    s.t.  0  <,   a  <,   C-max         in  a 

line  search,  yielding  am.  Let  w  <—  w  +  amd  and 

go  to  Step  1 . 
The  convergence  of  this  method  may  not  be  satisfactory 
since  it  represents  a  greedy  descent  method  in  the  reduced 
space  of  the  nonbasic  variables.  ReJclaitis  et  al .  [Ref.  30] 
suggest  a  convergence  acceleration  technique  by  using  a  more 
effective  unconstrained  search  method  as  long  as  the  basic 
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variables  do  not  change.  This  is  especially  relevant  in  our 
case  where  the  prior  iterate  frequently  takes  the  largest 
weight  and  remains  the  basic  variable.  Therefore  we  apply  a 
conjugate  gradient  method,  that  has  proven  its  effectiveness 
in  unconstrained  optimization  and  is  easily  implemented. 
The  modified  direction  for  the  first  iteration  and  any 
iteration  following  a  basis  change  becomes  : 

-  r!^    if  w.  >  0  or  r^   <    0 
dkNj  =  (62) 

0      otherwise, 
and  for  all  other  iterations  : 

k         k  "  r*  "2        k-i 

d  N  =  "  ^  +     ~  rk-l  ,,2    d/    •  <«> 

The  direction  for  the  basis  variable  is  the  same  as  in  the 
standard  reduced  gradient  method.  The  use  of  the  conjugate 
gradient  provides  significant  improvements.  In  a  typical 
master  problem  with  only  four  extreme  points,  the  reduced 
gradient  method  used  500  iterations  to  obtain  a  solution 
that  was  still  13%  worse  than  a  solution  obtained  with  the 
conjugate  gradient  method  in  28  iterations. 

Each  line  search  requires  frequent  evaluations  of  the 
objective  function  and  consumes  a  fair  amount  of  computation 
time.  We  can  improve  this  by  modifying  the  objective 
function  to  include  only  those  arcs  in  the  barrier  term  that 
are  potentially  able  to  violate  the  joint  capacity  con- 
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straints  at  the  current  iteration.  These  arcs  are  easily 
identified  and  recorded  in  a  set  Jy  containing  all  arcs  that 
ever  had  a  violation  in  any  of  the  extreme  points  generated 
in  the  subproblems .  The  arcs  j  4  Jy  cannot  violate  the 
joint  capacity  constraints  in  a  solution  to  the  master 
problem  that  is  a  convex  combination  of  the  extreme  points 
and  the  feasible  prior  iterate.  Thus  we  establish  the 
barrier  only  on  a  reduced  subset  of  the  joint  capacity 
constraints.  This  set  is  updated  at  each  solution  of  a 
subproblem  in  case  a  new  arc  experiences  a  capacity  viola- 
tion. This  procedure  amounts  to  a  modified  barrier  func- 
tion. However,  in  a  finite  number  of  iterations,  the  number 
of  arcs  in  Jy  will  take  a  fixed  value  IJyl  ^  |J|  and  no  more 
arcs  will  be  added  to  Jy.  From  this  point  on  we  are  back  to 
a  pure  barrier  function  as  derived  in  Chapter  2  and  conver- 
gence of  the  algorithm  is  preserved. 

C.   LOWER  AND  UPPER  BOUNDS 

Since  we  will  rarely  be  able  to  find  the  optimal 
solution  to  (P)  within  a  reasonable  number  of  iterations,  we 
need  to  establish  bounds  on  the  optimal  solution. 

Lower  bounds  can  be  derived  from  the  Lagrangian  dual 
problem 


(LR(u-j))     min  (c  -  u-.  A)  x  +  u-^b^  ,  which  provides  a  lower 
X€F 

bound  on  (P)  for  any  fixed  u-^  <,   0. 
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Using  the  dual  estimates   u-,  =  -  ji-,,  v.i.z. 

Uj  =  -  n  (b^  -  Ax)    for  each  solution  x  provided  by  the 

master  problem,  we  obtain 

(LR(|11))         min   (c  +  j^A)  x  -  (i1b1  .  (64) 

xeF 

Recalling  the  subproblem  (SP2k) :   min  (c  +  (La)  x  ,  we  find 

xeF 

by  comparison  that  both  objective  functions  differ  only  by 

the  constant  term  M^b-^  and  yield  the  same  optimal  solution 

k.  A  k- 

x  .   Thus,  we  obtain  a  valid  lower  bound  V(|i.,  *■)  by  subtract- 

ing  the  constant  term  M-ib-,  : 

V(G1k)  =  (c  +  G1kA)xk  -  M-1kb1.  (65) 

Furthermore,  since  in  the  limit  x   — »  x   and   |J,   — »  \L*    ,     it 

must  be  that  V((i1k)  — >    (c  +  ^   A)  x  -  |i1  b-,    =   ex  . 

Upper  bounds  on  the  optimal  solution  are  generated  in 

each  master  problem,  since  we  restrict  the  solution  to  an 

interior,  feasible  point.   Thus,  if  xk  =  Xkwk  solves  (MP2k) , 


the  upper  bound 


A],  Al, 

V(xK)  =   cxK  (66) 


is  readily  available  at  each  master  problem  solution.   Due 
to  the  convergence  property  (41)  of  the  barrier  function, 
it  must  be  that   V(x  )  — »  ex   as  h  — >   0. 

However,  we  are  usually  able  to  obtain  a  better  upper 
bound   by   utilizing   the   capacity   allocation   mechanism 

/NAT.  A], 

again.   Let  y  =  CA(x  )  at  some  iteration  k  and  let  x   be  an 
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interior  point.   Then,  for  all  arcs  j  e  Jv,   (Ax  -  b,)  •  <  0 

a       A    lc 

and      y  •  >  x  •    for  all  p, j   from  equations  (24).   Then 

H  J      r  J 

the  following  relationship  must  hold: 

V(x  )  =  ex     >    min  ex    >    min  ex   =  V(y), 

xeF  xeF 


A  k  A 

x^x  x<y 


A         A  1_ 

since  y  >  x 


Thus,  solving  SP1 (y)  yields  a  feasible  solution  with  value 
V(y)   <  V(xk)  . 


D.   THE  ALGORITHM  RSD (B) 

The  algorithm  RSD (B)  using  a  barrier  function  decomposi- 
tion can  now  be  presented.  An  initial  lower  bound  is 
obtained  by  solving  the  problem  without  the  joint  capacity 
constraints  (2).  If  this  solution  is  feasible  in  (P),  it  is 
optimal.  Otherwise  we  obtain  a  feasible  solution  via  the 
capacity  allocation  mechanism,  which  gives  an  initial  upper 
bound.  The  algorithm  generates  a  sequence  of  extreme  points 
in  the  subproblem  and  an  interior  point  solution  in  the 
master  problem  until  €-optimality  is  achieved.  As  a 
heuristic,  we  invoke  the  capacity  allocation  mechanism  at 
every  r  iteration  of  the  master  problem  to  improve  the 
current  upper  bound  and  decrement  h  at  this  time,  even  if 
the  master  problem  is  not  solved  optimally  for  the  cur- 
rent h. 
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Notation  : 

X      matrix  of  retained  extreme  points  at  iteration  k 

Ak 

x     current  master  problem  solution 

xM     previous  master  problem  solution 

Ak         •    k 

X      matrix  X   augmented  with  xM 

k  k- 

x      optimal  solution  to  subproblem  (SP2  ) 

H(X)  convex  hull  of  X 

CA(x)  a  capacity  allocation  based  on  x 

Jy  set  of  joint  capacity  arcs  which  are  violated  in 

at  least  one  subproblem  solution  x  . 

h  barrier  parameter  used  in  1    parameter  update 

r  maximal  number  of  retained  extreme  points 

e  stopping  criteria  for  near-optimality . 

Algorithm  RSD (B) 

Input   :   The  network  G  =  {I,  J},  joint  capacity  vector  b-,, 
cost  vector  c  and  supply/demand  vector  b2 • 

Ak  — 

Output  :   Best  obtained  solution  x   and  final  bounds  V,  V. 

Step  0  :  (Initialization) 

Select  h°  >  0,  €  >  0,  r  >  1,  0  <  b  <  1. 

Set  k  =  0,  1=0, 

X°  =  0  ,  Jy  =  0  ,  ji-,9  =  0   for  all  j  e  J. 

Solve  (LR(0))  :  x°  =  argmin  {ex  |  x  €  F  }. 

Set  V  =  cx°.   Set  Jy  =  { j  |  (^  -Ax°) •  <  0}. 
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A  0 
If  Jv  =  0,  stop  with  x   optimal. 

Else,  set  y  =  CA(x°),  y'  =  b  *  y  and 

solve  SP1 (y' )  yielding  x  u  . 

Set  V  =  ex  °.   If  (V  -  V)/V  <  €,  exit  with  x  °,  V,  V, 

Else,  set  X1  =  0  ,  xM  =   x1  =  x  °  ,  k  -  1. 


Step  1  :  (Solve  SP2k) 

Solve  xk  =  argmin  {  (c  +  jl^AJx  |  x  e  F}, 
x 

where   (l^   =  h1  (b1  -  Axk)  "-1   for  all  j  e  Jv  and 

\lA   =   0  for  all  j  £    Jv. 

Set  Jv  =  Jv  U  {j  |  (b1    -   Axk)  <  0}. 

A    1_  A    1^      1-        A    },. 

Set  V(p.1K)  =  (c  +  \i1K)xK   -  \i±    b1. 

If  Vf^*)  <  V,  set  V  =  V(ii1lc)  . 

If  (V  -  V)  /  V  <  e,  exit  with  x*,  V,  V. 

If  (c  +  M-^A)  (xk  -  xk)  >  0,  xk  solves  B(h1,x)  . 

Set  h1  +  1  =  a  *  h1   (0  <  a  <,   1)  . 

Set  1=1+1  and  go  to  Step  2 . 
Else 

(i)   if  |Xk|  <  r  ,  Xk+1  =  Xk  U  {xk}. 

(n)  if  |X  |  =  r  ,  drop  the  column  of  X  which  had 

Ak  . 
the  smallest  weight  w£  in  the  convex  com- 

Al.  _  1_ 

bination  forming  x   and  replace  it  with  x  . 
Go  to  Step  2 . 
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Step  2  :  (Solve  the  master  problem  MP2k) 
Set  Xk+1  =  Xk+1  U  {xM} . 
Find   xk+1  =  argmin  {  B(h1fx)  |  x  €  H(Xk+1)  } 


=  Ii  w^+1xi  where   1  <  i  <  |Xk+1|  and 

.  t-h  Ak  +  1 

Xj_  is  the  1    column  in  X    . 
Set  xM  =  xk+1. 

If  cxK+i  <  V  ,  set  V  =  cxK   . 
If  (V  -  V)  /  V  <  e,  exit  with  x*1,  V,  V. 
If  k+1  is  an  integer  multiple  of  r,  do  a  capacity 

allocation  :  Set  y  =  CA(xk+1)  and  solve  SP1 (y) , 

yielding  V(y)  =  ex  k+1 . 

If  V(y)  <  V,  set  V  =  V(y) . 


If  (V  -  V)  /  V  <  e,  set  x3^  =  x  k+1  , 

exit  with  xk,  V,  V. 
Set  h1  +  1  =  a  *   h1   (0  <  a  <,   1)  and  1  =  1  +  1 
Set  k  =  k  +  1  and  go  to  Step  1 . 


C.   IMPLEMENTATION  OF  RSD (B) 

The  algorithm  RSD (B)  has  been  coded  in  FORTRAN  which  is 
still  an  extremely  efficient  language  for  mathematical 
programming  purposes  (see  MacLennan  [Ref.  31]) .  A  sophis- 
ticated data  structure  is  used  for  the  storage  of  sparse 
matrices  and  vectors.  Allowing  direct  communication  with 
the  X-system  solver  of  Brown  and  Graves  [Ref.  32].   Integer 
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arithmetic  is  performed  to  the  extent  possible,  especially 
for  the  subproblems  which  use  a  GNET  solver  [Ref.  1]  to 
produce  rapid  solutions. 

A  design  goal  was  set  to  provide  a  decomposition 
procedure  which  solves  only  a  single  commodity  problem  at  a 
time,  and  thus  operates  easily  within  a  modest  memory 
region  (say,  two  megabyte  virtual  storage  capacity) .  Other 
information  such  as  current  incumbent,  previous  iterate, 
prior  extreme  points,  etc.,  have  to  be  kept  on  external 
storage  devices.  This  approach  leads  to  considerable 
input/output  operations  at  the  expense  of  computation  time, 
but  the  maximum  problem  size  in  terms  of  number  of  com- 
modities remains  independent  of  the  available  virtual 
storage.  Also,  the  resulting  algorithm  is  highly  parallel 
by  commodities . 

The  implementation  of  the  master  problem  contains 
several  parameters  such  as  the  number  of  retained  extreme 
points,  stopping  criteria  for  optimality  at  each  iteration, 
the  final  interval  of  uncertainty  in  the  line  search  and  an 
upper  limit  on  the  maximum  number  of  line  searches  con- 
ducted. Furthermore,  the  weights  of  the  extreme  points  and 
the  objective  function  evaluation  are  subject  to  roundoff 
errors.  For  a  sensitive  objective  function,  special  care 
is  necessary  to  insure  convergence.  Fortunately,  we  will 
confirm  that  the  barrier  function  decomposition  is  a  robust 
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procedure  that  does  not  require  an  optimal  solution  to  the 
master  problem  at  each  iteration.  Good  results  are  obtained 
over  a  relative  broad  range  of  parameters. 

The  comparison  of  the  barrier  algorithm  RSD (B)  with  the 
penalty  algorithm  RSD(P)  of  Staniec  [Ref.  4]  reveals  that 
they  are  very  similar  in  their  sequential  structure.  This 
similarity  permits  embedding  both  algorithms  in  the  same 
computer  program  and  creates  the  potential  for  devising  a 
hybrid  algorithm  which  takes  advantage  of  each.  The 
relationship  between  RSD (B)  and  RSD(P)  is  easily  establish- 
ed. The  dual  estimates  obtained  from  the  penalty  approach 
take,  for  some  vector  x  and  joint  capacity  constraint  j ,  the 
form 

(L^  =  h  [(Ax  -  b1)^"1]j   ,  t  >  1  (67) 

versus 

£lj  =  h[(bl  "  Ax>-1]j  '  <bl  "  Ax>j  >  0  ,j  e  Jv    (68) 
for  the  barrier  approach.     However,  the  gradient  of  the 
penalty  function  at  some  vector  x  takes  the  same  form  as  for 
the  barrier  function,  namely 

v-Q(h,x)  =  c  +  (j|a  (69) 

versus 

VB(h,x)  =  c  +  ^A.  (70)  . 

Thus,  besides  the  fact  that  the  barrier  function  decomposi- 
tion requires  an  initial  interior  point,  both  methods  use 
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the  same  subproblems  and  master  problem  with  different  input 
for  the  dual  estimates  and  different  function  evaluation 
routines  in  the  line  search. 
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IV.  COMPUTATIONAL  RESULTS 

In  order  to  assess  the  capabilities  of  the  algorithm 
RSD(B),  we  solve  different  versions  of  the  test  problem 
described  in  Chapter  I.  This  problem  suite  has  been  exten- 
sively studied  by  Staniec  [Ref.  4],  using  different  al- 
gorithms. The  optimal  solution  for  four  and  ten  commodity 
problems  is  available  for  comparisons,  obtained  by  solving 
the  problem  (P)  using  the  X-system  [Ref.  32]  .  The  four 
commodity  problem  (4H)  has  approximately  13,200  constraints, 
41,600  variables  and  optimal  solution  value  130,739,585.  The 
ten  commodity  problem  (10H)  has  about  33,000  constraints, 
104,000  variables  and  optimal  solution  value  169,532,339. 
For  purposes  of  direct  comparisons,  the  penalty  algorithm 
RSD(P)  of  Staniec  [Ref.  4]  has  been  converted  into  our  data 
structure  and  uses  our  computer  program  framework.  RSD(P) 
has  been  improved  by  taking  advantage  of  the  conjugate 
gradient  modifications  in  the  master  problem. 

First,  we  will  evaluate  the  performance  of  RSD(B)  with 
different  initial  parameters  in  solving  problem  4H.  Subse- 
quently, the  comparisons  between  RSD(P)  and  RSD(B)  are 
presented  for  problems  4H  and  10H.  Possible  modifications  of 
algorithm  RSD(B)  are  then  discussed.  Finally,  we  test  with 
a  100  commodity  problem  having  more  than  one  million  vari- 
ables and  300,000  constraints. 
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A.  PERFORMANCE  OF  THE  ALGORITHM  RSD (B) 

The  following  results  were  obtained  on  an  IBM  3033  AP 
under  VM/CMS.   The  Central  Processor  Unit  (CPU)  utilization 
time  is  used  as  a  performance  measurement.  An  "e -optimality 
gap"  is  computed  from  the  current  upper  and  lower  bounds  as 
(  V  -  V  )  /  ex   or  estimated  as  (  V  -  V  )  /  V  . 

Since   the   barrier   function   decomposition   requires   an 
interior   starting   point,   we   will   first   investigate   the 
response  of  RSD (B)  to  different  starting  points.    Let  c 
denote  the  maximum  cost  value  in  the  network.   While  fixing 
h   =  2  *  c__v,  a  =  0.5,  and  r  =  7  in  problem  4H,  the  choice 

A  A 

of  the  parameter  b  establishing  y'  =  b  *  y  resulted  in  the 
optimality  gaps  listed  in  Table  1: 

Gap  (%) 


initial  gap 

66.76 

51.93 

39.42 

26.95 

iteration 

7 

10.21 

9.87 

9.56 

10.49 

iteration 

16 

4.60 

3.98 

4.60 

5.06 

iteration 

24 

3.35 

2.25 

3.05 

2.73 

iteration 

32 

2.03 

1.42 

1.53 

1.86 

iteration 

40 

1.19 

0.99 

1.28 

1.38 

iteration 

48 

0.75 

0.68 

0.89 

0.89 

reduction 

b 

0.80 

0.85 

0.90 

0.95 

Table  1  :  Response  to  Different  Starting  Points 
Problem  4H 
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If  the  parameter  b  is  selected  too  low,  the  dual  es- 
timates associated  with  the  resulting  interior  starting 
point  may  not  generate  good,  initial  extreme  points.  If  the 
starting  value  is  chosen  too  high,  the  interior  point  gets 
too  close  to  the  boundary  of  the  constraint  region.  The  sub- 
sequent line  search  in  the  master  problem  is  confined  to  a 
smaller  search  space  of  the  constraint  region  resulting  in 
slower  convergence.  Both  values,  b  =  0.85  and  b  =  0.9,  work 
well  in  the  test  problem.  Further  analysis  will  be  based  on 
b  =  0.9.  Good  results  with  RSD  were  obtained  for  any  value 
of  r  between  6  and  8. 

The  lower  bounds  obtained  with  the  algorithm  RSD (B)  are 

sensitive  to  the  initial  parameter  h  .   Recall  that  the  dual 

A  -1 

estimates  are  computed  as  ji-,  =  h(Ax-b-i)   .   If  the  barrier 

parameter  h  is  relatively  small,  we  approach  the  boundary 
very  soon.  As  the  slack  on  some  jointly  capacitated  arcs 
almost  vanishes,  its  reciprocal  may  take  huge  values  result- 
ing in  extreme  dual  estimates.  This  situation  leads  to  large 
oscillations  in  subsequent  solutions  of  the  subproblem.  The 
phenomenum  is  demonstrated  in  Figure  3  where  h  =  0.5  *  cmax* 
As  h  is  increased,  this  effect  seems  to  disappear,  as  shown 
in  Figures  4  through  6.  It  is  interesting  to  observe  that 
some  good  lower  bounds  are  obtained  under  all  conditions. 
The  extreme  result  obtained  for  the  small  value  of  h 
suggests  that  h°  >  cmax  is  the  better  choice. 
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Figure    3:    V(^)    for   hu   =   0.5*cmax,    Problem   4H 
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Figure  4:  V^)  for  h°  =  l*cmax,  Problem  4H 
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Figure    5:    V(fi1)    for   hu    =    l-5*cmax,    Problem   4H 
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Figure  6:  V(jj.1)  for  hu  =  2*cmax,  Problem  4H 
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If  we  want  to  select  an  initial  value  h  that  balances 
the  linear  part  ex  and  the  barrier  term  in  the  objective 
function,  we  have  to  take  into  account  that  I  Jy  I  is  not  fixed 
but  initially  increasing.  We  do  not  have  sufficient  inform- 
ation about  the  size  of  the  barrier  term  unless  at  least  one 
pilot  run  has  been  conducted.  Using  the  information  that 
about  6%  of  the  arcs  are  finally  contained  in  the  set  Jy  for 
problem  4H,  we  would  obtain  a  value  of  h  <  0.5  *  c  ,  which 
is  not  the  best  choice,  but  may  serve  as  a  lower  bound  on  h  . 

The  differences  obtained  for  the  upper  bounds  are  less 

Ak 

significant.   Figure   7   shows   the   sequence   of  values   ex 

obtained  for  the  same  choices  of  h   as  before. 
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Figure  7:  Response  of  cx^  to  h  ,  Problem  4H 
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Any  decrease  in  h  creates  a  disturbance  in  the  dual  es- 
timates. Since  the  master  problem  is  not  solved  optimally 
for  each  value  of  h,  the  algorithm  may  not  always  provide 
better  lower  bounds  before  the  next  update  of  h  occurs.  We 
are  still  able  to  get  good  lower  bounds  when  a  relatively 
moderate  rate  of  decrease  a  =  0.5  is  used. 

A  solution  trajectory  for  problem  4H  is  given  in  Figure  8 

using  h   =  1.5  *  c__v,  where  a  0.6%  optimal  solution  is 

in  d  4. ». 

obtained   within   200   seconds.     Instead   of   plotting   the 

current  upper  and  lower  bounds,  the  current  values  of  V(ji) 

Ak  . 

and  ex   are  shown.   We  observe  the  almost  strictly  decreasing 

Ak 

sequence  of  ex   although  the  master  problem  is  not  solved 

optimally.  The  improved  upper  bounds  obtained  by  the 
capacity  allocation  mechanism  are  denoted  as  ex'  and  indi- 
cate that  significant  gains  are  obtained  only  at  the  initial 

Ak 
stage.   The  differences  between  ex   and  ex'   nearly  vanish 

after  about   100   seconds   although  slightly  tighter  upper 

bounds  are  produced  throughout. 
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Figure  8:  Trajectory  of  RSD (B) ,  Problem  4H 

It  seems  to  be  typical  for  the  barrier  function  that  good 
upper  bounds  are  obtained  at  an  early  stage,  whereas  the 
lower  bounds  trail  behind. 

The  trajectory  of  the  objective  value  B(h,x  )  together 
with  its  linear  part  ex  is  given  in  Figure  9.  We  find  that 
B(h,x  )  still  approaches  the  optimum  ex  from  below,  the 
barrier  term  takes  negative  values  over  the  whole  range.  The 
steps  in  its  trajectory  are  due  to  the  decreases  of  h  by  a 
factor  of  0.5. 
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Figure  9:  Values  of  the  Objective  Function,  Problem  4H 

A  direct  comparison  between  the  algorithms  RSD(P)  and 
RSD  (B)  shows  that  the  barrier  function  decomposition  is  less 
effective  in  the  solution  of  the  smaller  problem  4H  but  is 
competitive  in  the  problem  10H.  The  solutions  to  problem  4H 
are  compared  in  Figure  10  were  h°  -  0.001*cmax  for  RSD(P)  and 
h°  =  l-5*cmax  for  RSD(B) . 
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Figure  10:  Comparison  of  RSD(P)  to  RSD (B) ,  Problem  4H 

If  we  use  the  same  initial  parameters  for  the  solution  of 
problem  10H,  we  obtain  an  interesting  result.  The  trajec- 
tories are  given  in  Figure  11.  Obviously,  the  initial 
barrier  and  penalty  parameters,  respectively,  are  too  low  in 
both  cases,  yielding  poor  and  oscillating  lower  bounds  for 
RSD(B)  versus  bad  upper  bounds  for  RSD(P)  due  to  insufficient 
penalty  on  the  capacity  violations.  (This  situation  suggests 
investigation  of  a  hybrid  algorithm,  incorporating  both 
RSD (B)  and  RSD(P) ) . 
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Figure  11:  RSD(P)  Versus  RSD (B)  ,  Problem  10Hf 
Same  Initial  Parameters 

The  initial  parameter  h  has  been  increased  until  good 
results  are  obtained  for  both  methods.  The  result  is  shown 
in  Figure  12.  We  observe  that  the  final  values  of  V(|i  ) 
decay  for  both  algorithms.  We  do  not  generate  improving 
lower  bounds.  Apparently,  the  parameter  h  is  updated  to 
rapidly  before  the  master  problem  generates  good  dual 
estimates.  Staniec  [Ref.  4]  reported  similiar  poor  lower 
bounding  for  problem  10H. 
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Figure  12:  RSD(P)  Versus  RSD(B),  Problem  10H, 
Best  Parameters 

In  summary,  we  find  that  RSD(B)  provides  good  upper 
bounds  at  each  iteration  and  does  not  depend  on  the  capacity 
allocation  routine  as  the  penalty  algorithm  does.  On  the 
other  hand,  the  penalty  algorithm  provides  better  initial 
dual  estimates,  resulting  in  better  lower  bounds.  Both 
algorithms  converge  to  a  good  solution  within  a  reasonable 
amount  of  time. 
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B.   POSSIBLE  MODIFICATIONS   OF  THE  ALGORITHM 

A  first  modification  could  be  to  postpone  a  further 
decrease  of  the  barrier  function  parameter  h  until  at  least 
one  better  lower  bound  has  been  obtained  compared  to  the 
bound  generated  for  the  previous  value  of  h.  Note,  however, 
that  the  sequence  of  upper  bounds  is  almost  strictly  decreas- 
ing from  iteration  to  iteration  and  controlled  by  the 
magnitude  of  h.  A  higher  barrier  parameter  results  in  higher 
upper  bounds.  Therefore,  this  choice  should  depend  on  the 
problem  at  hand.  If  a  fast  approximation  of  the  solution  is 
desired,  h  should  be  decreased  more  rapidly.  If  a  higher 
resolution  is  required  and  sufficient  computer  resources  are 
available,  a  better  solution  of  the  master  problem  is 
necessary,  yielding  better  intermediate  lower  bounds. 

Besides  such  a  change  of  input  parameters,  a  structural 
modification  was  also  investigated.  Since  each  capacity 
reallocation  routine  creates  a  subproblem  solution  that  is 
feasible  and  at  least  as  good  as  the  current  upper  bound, 
this  information  can  be  passed  to  the  master  problem  by  using 
this  solution  as  an  additional  "extreme"  point.  Experimenta- 
tion on  the  test  problem  showed  however  that  no  significant 
improvements  were  obtained.  This  may  be  because  the  dif- 
ference between  the  interior  point  solution  of  the  master 
problem   and   the   solution   of   the   capacity   reallocation 
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procedure  vanishes  while  the  lower  bounds  are  still  trailing 
behind.  Experimentation  on  other  problems  may  yield  dif- 
ferent results. 

The  capacity  allocation  mechanism  can  be  modified 
further.  For  RSD(B),  capacity  allocation  for  an  interior 
point  amounts  to  a  redistribution  of  the  available  slack.  As 
stated  in  equation  (24),  each  commodity  receives  the  same 
additional  amount.   A  proportional  allocation  would  be  given 

by       ypj  =  xpj  +  xpj  (bx  -  Ax)j  /  (a£)j.  (71) 

The  disadvantage  of  this  procedure  is  that  a  commodity  with 
zero  flow  on  that  arc  gets  zero  capacity  allocated.  To 
overcome  this  drawback,  a  convex  combination  of  both  methods 
is  possible: 

*PJ  =  *PJ  +  Bl  *PJ  (bl  "  A*>j  '    (A*>J 

+  A2  (b1    -  Ax) j  /  |P|  (72) 

where  fl-,,  fi»2  >  0  and  B-,  +  15^  =  !• 

Experimentation  with  the  test  problem  4H  showed  improve- 
ments in  the  upper  bounds  in  both  cases,  but  since  the  lower 
bounds  are  still  weak,  the  overall  gain  is  not  significant 
for  this  problem. 

RSD(B)  would  be  superior  to  RSD(P)  if  we  could  establish 
better  lower  bounds.  The  dual  estimates  (1  are  a  function  of 
the  barrier  parameter  h  and  the  slack  on  the  joint  capacity 
arcs.    Different  dual  estimates  are  obtained  if  we  use  a 
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different  value  of  h  for  each  individual  joint  capacity  arc. 
This  approach  leads  to  the  weighted  barrier  function  like  the 
one  proposed  by  Megiddo  [Ref.  14]  in  general  linear  programm- 
ing : 

B(h,x,w)  =  ex  -  h  I-  w •  ln(b1  -  Ax)  •  ,  (73) 

where  w  is  any  real  vector  with  positive  components.  Some 
experiments  were  done  with  different  weights.  All  arcs  that 
are  violated  in  the  initial  solution  with  completely  relaxed 
joint  capacity  constraints  are  more  likely  to  be  tight  in  the 
optimal  solution.  Therefore  they  are  assigned  a  reduced 
weight  to  enable  a  faster  approach  to  the  boundary.  Another 
weight  factor  used  was  proportional  to  the  remaining  slack  on 
the  corresponding  arc.  Unfortunately,  neither  attempt 
provided  better  solutions. 

The  final  modification  is  a  hybrid  barrier-penalty 
algorithm.  Starting  with  the  barrier  function  decomposition 
in  problem  10H,  we  shift  to  the  penalty  algorithm  after  16 
iterations  using  the  extreme  points  generated  so  far.  The 
main  problem  is  the  adjustment  of  the  parameter  h.  Since  we 
obtain  good  lower  bounds  for  relatively  small  values  of  h, 
we  reset  h  to  the  same  initial  value  that  we  used  in  the 
independent  approach.  The  results  displayed  in  Figure  13  are 
unimpressive  but  further  experimentation  may  improve  results. 
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Figure  13  :  Trajectory  of  the  Hybrid  Method 

Finally,  we  obtain  a  solution  to  the  100  commodity 
network  flow  problem  100H  which  is  about  3.5%  optimal  within 
2650  seconds  after  25  iterations  and  about  2.5%  optimal 
within  3800  seconds  after  38  iterations.  The  initial  value 
h  has  been  increased  to  h  =  ^*cmax*  Staniec  [Ref.  4] 
reports  a  solution  to  this  problem  obtained  by  a  penalty 
method  that  achieved  4  %  optimality  after  1000  seconds  and 
finished  with  1.5%  optimality  after  3000  seconds.  Our  method 
did  not  show  the  same  performance.  The  solution  trajectory 
as  shown  in  Figure  14  seems  to  indicate  that  the  upper  bounds 
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are  already  very  tight  whereas  the  lower  bounds  are  again 
poor.  The  objective  function  value  (not  shown)  is  still  less 
than  the  lower  bound,  a  further  decrease  of  h  is  necesarry 
for  convergence. 
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Figure  14  :  Solution  Trajectory,  Problem  100H 

The  distribution  of  the  CPU  time  between  the  master 
problem  and  the  subproblems  changes  as  we  increase  the  number 
of  commodities.  Although  the  solution  of  each  pure  network 
flow  problem  requires  less  than  one  second  CPU  time  per 
commodity  in  the  average  for  all  test  problems,  we  find  that 
the   subproblems   are  expensive  to   solve  and  consume  the 
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largest  portion  of  the  total  CPU  time.  This  distribution  is 
shown  in  Table  2.  The  amount  of  time  not  used  by  the  master 
problem  or  the  subproblem  is  mostly  consumed  in  the 
input/output  operations. 


Problem  4H    Problem  10H    Problem  100H 


Master  Problem 

25.5% 

11.2% 

7.1% 

Subproblem 

49.9% 

66.9% 

86.3% 

Other 

24.6% 

21.9% 

6.6% 

Table  2:   Distribution  of  CPU  Time  for  Test 
Problems  of  Increasing  Size 

However,  the  subproblem  solves  for  each  commodity  indepen- 
dently and  more  than  one  pure  network  flow  problem  could  be 
solved  at  a  time.  This  feature  makes  our  method  highly 
parallel  with  a  potential  reduction  of  nearly  1/|P|  in 
elapsed  computation  time  for  the  subproblem. 
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V.  CONCLUSIONS 

The  method  of  logarithmic  barrier  function  decomposition 
is  a  useful  tool  for  the  solution  of  large-scale  multicom- 
modity  flow  problems.  The  algorithm  RSD (B)  is  a  variant  of 
the  price-directive  decomposition  method.  It  generates  a 
sequence  of  interior  points  which  provide  intermediate 
feasible  solutions  while  converging  towards  the  optimum.  The 
technique  of  restricted  simplicial  decomposition  (RSD)  proves 
to  be  useful  in  the  solution  of  the  nonlinear  master  problem. 
However,  RSD(B)  seems  to  be  robust  and  does  not  require  an 
optimal  solution  to  the  master  problem.  It  is  competitive 
with  penalty  decomposition  methods  and  relatively  easy  to 
implement.  Since  it  achieves  good  feasible  solutions  at  an 
early  stage,  it  may  be  used  as  a  starting  technique  in  other 
algorithms  like  the  hybrid  barrier-penalty  technique.  The 
use  of  RSD(B)  in  other  large-scale  multicommodity  flow 
problems  is  recommended. 
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